home *** CD-ROM | disk | FTP | other *** search
/ Cream of the Crop 26 / Cream of the Crop 26.iso / os2 / octa209s.zip / octave-2.09 / libcruft / blas / dtrmm.f < prev    next >
Text File  |  1997-01-29  |  11KB  |  356 lines
  1.       SUBROUTINE DTRMM ( SIDE, UPLO, TRANSA, DIAG, M, N, ALPHA, A, LDA,
  2.      $                   B, LDB )
  3. *     .. Scalar Arguments ..
  4.       CHARACTER*1        SIDE, UPLO, TRANSA, DIAG
  5.       INTEGER            M, N, LDA, LDB
  6.       DOUBLE PRECISION   ALPHA
  7. *     .. Array Arguments ..
  8.       DOUBLE PRECISION   A( LDA, * ), B( LDB, * )
  9. *     ..
  10. *
  11. *  Purpose
  12. *  =======
  13. *
  14. *  DTRMM  performs one of the matrix-matrix operations
  15. *
  16. *     B := alpha*op( A )*B,   or   B := alpha*B*op( A ),
  17. *
  18. *  where  alpha  is a scalar,  B  is an m by n matrix,  A  is a unit, or
  19. *  non-unit,  upper or lower triangular matrix  and  op( A )  is one  of
  20. *
  21. *     op( A ) = A   or   op( A ) = A'.
  22. *
  23. *  Parameters
  24. *  ==========
  25. *
  26. *  SIDE   - CHARACTER*1.
  27. *           On entry,  SIDE specifies whether  op( A ) multiplies B from
  28. *           the left or right as follows:
  29. *
  30. *              SIDE = 'L' or 'l'   B := alpha*op( A )*B.
  31. *
  32. *              SIDE = 'R' or 'r'   B := alpha*B*op( A ).
  33. *
  34. *           Unchanged on exit.
  35. *
  36. *  UPLO   - CHARACTER*1.
  37. *           On entry, UPLO specifies whether the matrix A is an upper or
  38. *           lower triangular matrix as follows:
  39. *
  40. *              UPLO = 'U' or 'u'   A is an upper triangular matrix.
  41. *
  42. *              UPLO = 'L' or 'l'   A is a lower triangular matrix.
  43. *
  44. *           Unchanged on exit.
  45. *
  46. *  TRANSA - CHARACTER*1.
  47. *           On entry, TRANSA specifies the form of op( A ) to be used in
  48. *           the matrix multiplication as follows:
  49. *
  50. *              TRANSA = 'N' or 'n'   op( A ) = A.
  51. *
  52. *              TRANSA = 'T' or 't'   op( A ) = A'.
  53. *
  54. *              TRANSA = 'C' or 'c'   op( A ) = A'.
  55. *
  56. *           Unchanged on exit.
  57. *
  58. *  DIAG   - CHARACTER*1.
  59. *           On entry, DIAG specifies whether or not A is unit triangular
  60. *           as follows:
  61. *
  62. *              DIAG = 'U' or 'u'   A is assumed to be unit triangular.
  63. *
  64. *              DIAG = 'N' or 'n'   A is not assumed to be unit
  65. *                                  triangular.
  66. *
  67. *           Unchanged on exit.
  68. *
  69. *  M      - INTEGER.
  70. *           On entry, M specifies the number of rows of B. M must be at
  71. *           least zero.
  72. *           Unchanged on exit.
  73. *
  74. *  N      - INTEGER.
  75. *           On entry, N specifies the number of columns of B.  N must be
  76. *           at least zero.
  77. *           Unchanged on exit.
  78. *
  79. *  ALPHA  - DOUBLE PRECISION.
  80. *           On entry,  ALPHA specifies the scalar  alpha. When  alpha is
  81. *           zero then  A is not referenced and  B need not be set before
  82. *           entry.
  83. *           Unchanged on exit.
  84. *
  85. *  A      - DOUBLE PRECISION array of DIMENSION ( LDA, k ), where k is m
  86. *           when  SIDE = 'L' or 'l'  and is  n  when  SIDE = 'R' or 'r'.
  87. *           Before entry  with  UPLO = 'U' or 'u',  the  leading  k by k
  88. *           upper triangular part of the array  A must contain the upper
  89. *           triangular matrix  and the strictly lower triangular part of
  90. *           A is not referenced.
  91. *           Before entry  with  UPLO = 'L' or 'l',  the  leading  k by k
  92. *           lower triangular part of the array  A must contain the lower
  93. *           triangular matrix  and the strictly upper triangular part of
  94. *           A is not referenced.
  95. *           Note that when  DIAG = 'U' or 'u',  the diagonal elements of
  96. *           A  are not referenced either,  but are assumed to be  unity.
  97. *           Unchanged on exit.
  98. *
  99. *  LDA    - INTEGER.
  100. *           On entry, LDA specifies the first dimension of A as declared
  101. *           in the calling (sub) program.  When  SIDE = 'L' or 'l'  then
  102. *           LDA  must be at least  max( 1, m ),  when  SIDE = 'R' or 'r'
  103. *           then LDA must be at least max( 1, n ).
  104. *           Unchanged on exit.
  105. *
  106. *  B      - DOUBLE PRECISION array of DIMENSION ( LDB, n ).
  107. *           Before entry,  the leading  m by n part of the array  B must
  108. *           contain the matrix  B,  and  on exit  is overwritten  by the
  109. *           transformed matrix.
  110. *
  111. *  LDB    - INTEGER.
  112. *           On entry, LDB specifies the first dimension of B as declared
  113. *           in  the  calling  (sub)  program.   LDB  must  be  at  least
  114. *           max( 1, m ).
  115. *           Unchanged on exit.
  116. *
  117. *
  118. *  Level 3 Blas routine.
  119. *
  120. *  -- Written on 8-February-1989.
  121. *     Jack Dongarra, Argonne National Laboratory.
  122. *     Iain Duff, AERE Harwell.
  123. *     Jeremy Du Croz, Numerical Algorithms Group Ltd.
  124. *     Sven Hammarling, Numerical Algorithms Group Ltd.
  125. *
  126. *
  127. *     .. External Functions ..
  128.       LOGICAL            LSAME
  129.       EXTERNAL           LSAME
  130. *     .. External Subroutines ..
  131.       EXTERNAL           XERBLA
  132. *     .. Intrinsic Functions ..
  133.       INTRINSIC          MAX
  134. *     .. Local Scalars ..
  135.       LOGICAL            LSIDE, NOUNIT, UPPER
  136.       INTEGER            I, INFO, J, K, NROWA
  137.       DOUBLE PRECISION   TEMP
  138. *     .. Parameters ..
  139.       DOUBLE PRECISION   ONE         , ZERO
  140.       PARAMETER        ( ONE = 1.0D+0, ZERO = 0.0D+0 )
  141. *     ..
  142. *     .. Executable Statements ..
  143. *
  144. *     Test the input parameters.
  145. *
  146.       LSIDE  = LSAME( SIDE  , 'L' )
  147.       IF( LSIDE )THEN
  148.          NROWA = M
  149.       ELSE
  150.          NROWA = N
  151.       END IF
  152.       NOUNIT = LSAME( DIAG  , 'N' )
  153.       UPPER  = LSAME( UPLO  , 'U' )
  154. *
  155.       INFO   = 0
  156.       IF(      ( .NOT.LSIDE                ).AND.
  157.      $         ( .NOT.LSAME( SIDE  , 'R' ) )      )THEN
  158.          INFO = 1
  159.       ELSE IF( ( .NOT.UPPER                ).AND.
  160.      $         ( .NOT.LSAME( UPLO  , 'L' ) )      )THEN
  161.          INFO = 2
  162.       ELSE IF( ( .NOT.LSAME( TRANSA, 'N' ) ).AND.
  163.      $         ( .NOT.LSAME( TRANSA, 'T' ) ).AND.
  164.      $         ( .NOT.LSAME( TRANSA, 'C' ) )      )THEN
  165.          INFO = 3
  166.       ELSE IF( ( .NOT.LSAME( DIAG  , 'U' ) ).AND.
  167.      $         ( .NOT.LSAME( DIAG  , 'N' ) )      )THEN
  168.          INFO = 4
  169.       ELSE IF( M  .LT.0               )THEN
  170.          INFO = 5
  171.       ELSE IF( N  .LT.0               )THEN
  172.          INFO = 6
  173.       ELSE IF( LDA.LT.MAX( 1, NROWA ) )THEN
  174.          INFO = 9
  175.       ELSE IF( LDB.LT.MAX( 1, M     ) )THEN
  176.          INFO = 11
  177.       END IF
  178.       IF( INFO.NE.0 )THEN
  179.          CALL XERBLA( 'DTRMM ', INFO )
  180.          RETURN
  181.       END IF
  182. *
  183. *     Quick return if possible.
  184. *
  185.       IF( N.EQ.0 )
  186.      $   RETURN
  187. *
  188. *     And when  alpha.eq.zero.
  189. *
  190.       IF( ALPHA.EQ.ZERO )THEN
  191.          DO 20, J = 1, N
  192.             DO 10, I = 1, M
  193.                B( I, J ) = ZERO
  194.    10       CONTINUE
  195.    20    CONTINUE
  196.          RETURN
  197.       END IF
  198. *
  199. *     Start the operations.
  200. *
  201.       IF( LSIDE )THEN
  202.          IF( LSAME( TRANSA, 'N' ) )THEN
  203. *
  204. *           Form  B := alpha*A*B.
  205. *
  206.             IF( UPPER )THEN
  207.                DO 50, J = 1, N
  208.                   DO 40, K = 1, M
  209.                      IF( B( K, J ).NE.ZERO )THEN
  210.                         TEMP = ALPHA*B( K, J )
  211.                         DO 30, I = 1, K - 1
  212.                            B( I, J ) = B( I, J ) + TEMP*A( I, K )
  213.    30                   CONTINUE
  214.                         IF( NOUNIT )
  215.      $                     TEMP = TEMP*A( K, K )
  216.                         B( K, J ) = TEMP
  217.                      END IF
  218.    40             CONTINUE
  219.    50          CONTINUE
  220.             ELSE
  221.                DO 80, J = 1, N
  222.                   DO 70 K = M, 1, -1
  223.                      IF( B( K, J ).NE.ZERO )THEN
  224.                         TEMP      = ALPHA*B( K, J )
  225.                         B( K, J ) = TEMP
  226.                         IF( NOUNIT )
  227.      $                     B( K, J ) = B( K, J )*A( K, K )
  228.                         DO 60, I = K + 1, M
  229.                            B( I, J ) = B( I, J ) + TEMP*A( I, K )
  230.    60                   CONTINUE
  231.                      END IF
  232.    70             CONTINUE
  233.    80          CONTINUE
  234.             END IF
  235.          ELSE
  236. *
  237. *           Form  B := alpha*A'*B.
  238. *
  239.             IF( UPPER )THEN
  240.                DO 110, J = 1, N
  241.                   DO 100, I = M, 1, -1
  242.                      TEMP = B( I, J )
  243.                      IF( NOUNIT )
  244.      $                  TEMP = TEMP*A( I, I )
  245.                      DO 90, K = 1, I - 1
  246.                         TEMP = TEMP + A( K, I )*B( K, J )
  247.    90                CONTINUE
  248.                      B( I, J ) = ALPHA*TEMP
  249.   100             CONTINUE
  250.   110          CONTINUE
  251.             ELSE
  252.                DO 140, J = 1, N
  253.                   DO 130, I = 1, M
  254.                      TEMP = B( I, J )
  255.                      IF( NOUNIT )
  256.      $                  TEMP = TEMP*A( I, I )
  257.                      DO 120, K = I + 1, M
  258.                         TEMP = TEMP + A( K, I )*B( K, J )
  259.   120                CONTINUE
  260.                      B( I, J ) = ALPHA*TEMP
  261.   130             CONTINUE
  262.   140          CONTINUE
  263.             END IF
  264.          END IF
  265.       ELSE
  266.          IF( LSAME( TRANSA, 'N' ) )THEN
  267. *
  268. *           Form  B := alpha*B*A.
  269. *
  270.             IF( UPPER )THEN
  271.                DO 180, J = N, 1, -1
  272.                   TEMP = ALPHA
  273.                   IF( NOUNIT )
  274.      $               TEMP = TEMP*A( J, J )
  275.                   DO 150, I = 1, M
  276.                      B( I, J ) = TEMP*B( I, J )
  277.   150             CONTINUE
  278.                   DO 170, K = 1, J - 1
  279.                      IF( A( K, J ).NE.ZERO )THEN
  280.                         TEMP = ALPHA*A( K, J )
  281.                         DO 160, I = 1, M
  282.                            B( I, J ) = B( I, J ) + TEMP*B( I, K )
  283.   160                   CONTINUE
  284.                      END IF
  285.   170             CONTINUE
  286.   180          CONTINUE
  287.             ELSE
  288.                DO 220, J = 1, N
  289.                   TEMP = ALPHA
  290.                   IF( NOUNIT )
  291.      $               TEMP = TEMP*A( J, J )
  292.                   DO 190, I = 1, M
  293.                      B( I, J ) = TEMP*B( I, J )
  294.   190             CONTINUE
  295.                   DO 210, K = J + 1, N
  296.                      IF( A( K, J ).NE.ZERO )THEN
  297.                         TEMP = ALPHA*A( K, J )
  298.                         DO 200, I = 1, M
  299.                            B( I, J ) = B( I, J ) + TEMP*B( I, K )
  300.   200                   CONTINUE
  301.                      END IF
  302.   210             CONTINUE
  303.   220          CONTINUE
  304.             END IF
  305.          ELSE
  306. *
  307. *           Form  B := alpha*B*A'.
  308. *
  309.             IF( UPPER )THEN
  310.                DO 260, K = 1, N
  311.                   DO 240, J = 1, K - 1
  312.                      IF( A( J, K ).NE.ZERO )THEN
  313.                         TEMP = ALPHA*A( J, K )
  314.                         DO 230, I = 1, M
  315.                            B( I, J ) = B( I, J ) + TEMP*B( I, K )
  316.   230                   CONTINUE
  317.                      END IF
  318.   240             CONTINUE
  319.                   TEMP = ALPHA
  320.                   IF( NOUNIT )
  321.      $               TEMP = TEMP*A( K, K )
  322.                   IF( TEMP.NE.ONE )THEN
  323.                      DO 250, I = 1, M
  324.                         B( I, K ) = TEMP*B( I, K )
  325.   250                CONTINUE
  326.                   END IF
  327.   260          CONTINUE
  328.             ELSE
  329.                DO 300, K = N, 1, -1
  330.                   DO 280, J = K + 1, N
  331.                      IF( A( J, K ).NE.ZERO )THEN
  332.                         TEMP = ALPHA*A( J, K )
  333.                         DO 270, I = 1, M
  334.                            B( I, J ) = B( I, J ) + TEMP*B( I, K )
  335.   270                   CONTINUE
  336.                      END IF
  337.   280             CONTINUE
  338.                   TEMP = ALPHA
  339.                   IF( NOUNIT )
  340.      $               TEMP = TEMP*A( K, K )
  341.                   IF( TEMP.NE.ONE )THEN
  342.                      DO 290, I = 1, M
  343.                         B( I, K ) = TEMP*B( I, K )
  344.   290                CONTINUE
  345.                   END IF
  346.   300          CONTINUE
  347.             END IF
  348.          END IF
  349.       END IF
  350. *
  351.       RETURN
  352. *
  353. *     End of DTRMM .
  354. *
  355.       END
  356.